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We present a scheme for robust finite temperature quantum simulation of stabilizer Hamiltonians. 
The scheme is designed for realization in a physical system consisting of a finite set of neutral atoms 
trapped in an addressable optical lattice that are controllable via 1- and 2-body operations together 
with dissipative 1-body operations such as optical pumping. We show that these minimal physical 
constraints suffice for design of a quantum simulation scheme for any stabilizer Hamiltonian at 
either finite or zero temperature. We demonstrate the approach with application to the abelian and 
non-abelian toric codes. 



I. INTRODUCTION 



There has been much recent work on approaches to experimentally engineer many-body quantum phases of mat- 
ter [IHZ] ■ In particular there is a wide array of lattice Hamiltonians whose ground states are novel quantum phases that 
are not yet known to exist in natural systems [SHE]; to physically realize such phases, one must generate the models 
artificially. Additionally, many of these novel phases have potential applications to quantum information I13H15) . 
A robust experimental realization of such phases of matter might be a route to building a fault-tolerant quantum 
computer [TS] . 

One such class of lattice Hamiltonians are quantum stabilizer Hamiltonians [HI [T3l [T7H24] . Stabilizer Hamiltonians 
are composed of commuting multi-qubit Pauli operators and play a key role in quantum error correction. Encod- 
ing quantum information into ground states of stabilizer Hamiltonians provides a natural physical architecture for 
realization of quantum error correcting codes and quantum memory in qubit arrays (13) . Because of their character- 
istic degenerate ground states and finite energy gaps to local excitations, they offer physical protection for encoded 
quantum information. Much recent interest has therefore focused on the quantum phases generated by stabilizer 
Hamiltonians, their ability to generate topological order and their relation to novel exotic phases [71 |2"51 - |2"T| . 

While the quantum codes derived from ground states of stabilizer Hamiltonians are of broad theoretical interest in 
quantum information theory, the transition from theoretical characterization to actual realization of these in physical 
systems is nevertheless impeded by a number of difficulties. Many stabilizer Hamiltonians are composed of non-local 
or many-body terms that are difficult to realize in practice. This usually requires significant effort and associated 
overhead in quantum engineering of interactions. Both trapping of ions in Paul traps and of atoms or molecules in 
optical lattices offer the ability to generate the required many-body interactions. Additionally, a robust simulation 
of a lattice Hamiltonian requires an entropy sink to remove entropy that accumulates from noise in the control 
operations and interactions with the environment [28 30 . Entropy can be removed either actively via algorithmic 
error correction, or passively via an effective coupling to an external reservoir 26, 29, 31 33 . In this work we will 
show how to overcome these challenges for the robust quantum simulation of both abelian and non-abelian toric code 
Hamiltonians and generate an effective coupling to a low temperature thermal reservoir. 

In the rest of this paper we first provide an introduction to stabilizer Hamiltonians in Section [H] In Section III 



then summarize the properties of the well known abelian and non-abelian toric code Hamiltonians of Kitaev, which 
constitute key examples of stabilizer Hamiltonians with topologically ordered ground states and abelian and non- 
abelian excitations, respectively. We then present a detailed discussion of thermalization for stabilizer Hamiltonians 
in Section IV The analysis given here provides both a generalization and a more efficient approach to generating 
thermalization than that employed in our previous work [30j . This improved approach to thermalization constitutes 
the main new set of results in this work. Following this, we summarize several routes to simulation of the stabilizer 
Hamiltonian and components needed for thermalization, including a non-perturbative approach that is considerably 
more efficient than our previous perturbative stroboscopic approach (Section [Vb. We outline how such a thermal 
quantum simulation of the abelian toric code may be implemented with a finite set of neutral atoms trapped in 
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an optical lattice and indicate what is required in order to generalize this to thermal quantum simulation of the 



non-abelian toric code. Section VI concludes with a brief summary. 



II. STABILIZER HAMILTONIANS 

Consider a class of lattice systems whose degrees of freedom are d— level systems (qudits) . We consider Hamiltonians 
with local n-body interactions 

^ = -E^E^ a) 

where is a spatial neighborhood, v labels the type of interaction, and {J v } are the interaction strengths with J„ > 0. 
We assume that {hf^} are local n-body projection operators with maximum eigenvalue 1. We will now consider the 
eigenstates of a local term h 6 {hj^}. We will label the eigenstates of h as | e^} , where e > 0, the eigenvalue of h is 
(1 — e), and the label d distinguishes degeneracies: 

h\e d ) = (l-e)\e d ). (2) 



With this notation, the ground state of each local term in the Hamiltonian is |0), with eigenvalue 1 and e = 0. We 



can consider the eigenoperators {b d , b d } that span the eigenspace of h\ 



b\ \e c ) = S e>0 \e d ) , b d \e c ) = 6 d , c |0) (3) 

h = l-^2b% (4) 

d 

An arbitrary eigenoperator can be formed from products of {b d ,b d }. 

We will now consider the case of such Hamiltonians for which all the commute. These constitute the class 
of stabilizer Hamiltonians, whose ground states include the well known stabilizer codes (T3J [TTJ [53]. In this case, 
eigenstates of H will be simultaneous eigenstates of all hKr. In particular the ground states of H will be the ground 
state of all local terms: 

Jfr|*o> = |*o) V(.A/» (5) 
Excited eigenstates of H can be generate by applying products of b\ to the ground state. For example, the state 

\^e^,^d),(M\^d)]} = bl^ d bl t ^ d \^ ) (6) 

has purely localized quasiparticle excitations at neighborhoods Af and A/"'. Therefore, the eigenoperators of H are 
consequently purely local, and can be formed from local superpositions of products of the {b\f d }. This locality of the 
eigenoperators is a key feature of stabilizer-like Hamiltonians. We expect that an arbitrary translationally invariant 
local Hamiltonian with non-communting terms will have momentum eigenstates, and therefore the eigenoperators will 
be extensive superpositions of local operators. 

The energy cost of a local quasiparticle excitation (A/ - , fi,a) is AE% = e^J^. An arbitrary one-qudit operator 
acting on qudit j, o~j, will be formed from a finite local sum of products of the eigenoperators of all neighborhoods 
Af{,i € Mi- Therefore an arbitrary one-qudit operator will only create, annihilate or translate quasiparticles within 
the local region of neighborhoods Mi. In contrast, an arbitrary Hamiltonian with non-commuting local terms will have 
eigenstates with propagating quasiparticles; consequently a local operator acting on an eigenstates will generically 
create an non-local superposition of quasiparticles. 

III. TOPOLOGICAL PHASES AND TORIC CODES 

Topologically ordered phases of matter are 2D quantum liquid states with no broken conventional symmetry [151134) . 
These topological phases have a quantum ordering which cannot be detected by a local order parameter. The 
topological order results in a robust ground state degeneracy on surfaces with a non-trivial topology and a finite 
energy gap to anyonic quasiparticle excitations. Different ground states can only be distinguished by non-local 
operators that wind around a non-contractible loop of the surface. 
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Such topological phases have been proposed as the basis for a physically fault tolerant quantum computer [TTJ IT3T - 
fl5] . Quantum information can be encoded in the degenerate ground states; since all local operators cause transitions to 
excited state above the finite gap, these phases are relatively insensitive to local perturbations. Furthermore, both the 
tunneling amplitude between ground states and splitting of the ground state degeneracy are suppressed exponentially 
in the system size. Logical operations can be performed by creating and braiding the anyonic excitations [35] . 
Since the result of braiding operations depends only on the topological properties of the braids, not the precise 
details of the braiding paths, these logical control operations are intrinsically robust against noisy control operations. 
For certain phases with a sufficiently rich topological order, braiding operations form a universal set of quantum 

gates HD US ESI EH- 

To robustly physically realize such a robust topological phase, it is essential to maintain equilibrium with a low 
temperature external reservoir that can remove entropy and any associated accumulation of excitations due to envi- 
ronmental noise and/or noisy control operations [TH Q21 HH [2H1 13H1 ISH]- A generic quantum simulation (e.g., with 
trapped cold atoms) is an open quantum system that is not intrinsically in equilibrium with an external thermal 
reservoir. Thus maintaining thermal equilibrium at low temperatures requires explicitly generating an effective cou- 
pling to an external reservoir. It is important to note there that while theoretical studies have shown that in the 
thermodynamic limit, topological order is destroyed at any finite temperature [40] . in a finite sized system there is a 
finite temperature crossover below which the topological order is preserved [41]. 

Kitaev has introduced a class of exactly soluble lattice models with topologically ordered ground states [TTJ [22] . 
The Hilbert space of these systems is defined in general by a set of qudits that sit on the links of an oriented square 
lattice. When this model is placed on a lattice on a torus, this is known as the toric code Hamiltonian. We outline 
here the generic model of the abelian toric code as well as its non-abelian generalizations. Each qudit state is labeled 
by an element of a finite group G, such that the local Hilbert space on each link is {\g),g E G}. For each qudit, we 
define the operators that perform left and right multiplication 

L%\g) = \hg) L h _\g) = \gh- 1 ), (7) 

as well as the projection operators 

I$\g) = S h , g \g) T t l\g) = 5 h ^ g \g). (8) 
The toric code Hamiltonian is a stabilizer-like Hamiltonian comprising commuting 4-body interactions: 

H G = -Y J A V ~Y J B P , (9) 

v p 

where the A v 's are 4-body interactions defined on the vertices of the lattice and the B p 's are 4-body interactions 
defined on the plaquettes of the lattice. The vertex terms are given by 

^41^^-/-,. ( io ) 

1 1 gee 

where {i, j, k,l} € v and are ordered clockwise around v. The vertex operators can be viewed as gauge transformations, 
and eigenvalue 1 eigenstate of A v are considered gauge invariant. Any state for which A„|\I/o) l^o) on a given vertex 
is therefore not gauge invariant, and considered to have a non-trivial electric charge at vertex v. 
The plaquette operators are given by: 

B p = Yl T a J.T a J.T a + \T a + \, (11) 

SlS2S3S4 = l 

where the sum over {31,32,53,54} whose product is the identity. Any state which is not a 1 eigenvalue eigenstate of 
some Bp is considered to have a non-trivial magnetic charge on the face of plaquette p. 

All the A v and B p commute, and so the ground state of Hq is a simultaneous eigenstate of all A v and B p operators 
with eigenvalue 1: 

A v |*o) = |*o) , B p |¥ ) = |*o) ■ (12) 

The ground state has no electric and magnetic charges, and there is a finite gap to electric and magnetic quasiparticle 
excitations. The stabilizer-like form of H means that all excited states have purely localized electric or magnetically 
charged quasiparticle excitations. The spectrum of H includes electric charges on the vertices, magnetic charges on 
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the plaquettes, and bound dyonic state of electric and magnetic charges. 

Magnetic charges are labeled by the conjugacy classes of G, where the conjugacy class Ch of element h is defined as 

Ch^ighg-^geG}. (13) 

The centralizers of each conjugacy class Ch, which commute with all elements of Ch label the electric charges. Pairs 
of neighboring excitations can be created by applying one-qudit operators to a link. For example, to create a pair of 
magnetic fluxes on neighboring plaquettes, we can define an operator E pp ,([h]) which acts on the link connecting the 
plaquettes p and // [27] : 

m^mg 1 ) = E+ p ,([h}) |* ) , E+ p ,([h}) = -±= £ L h _, (14) 

V \i n i \ h£[h] 

where \ra p h \rr^ p h )) is a state with magnetic fluxes of charge [h] on plaquettes p,p'. The operator E pp ,([h]) is a one 
qudit operator acting on the link connecting p and p' . An arbitrary one-qudit operator applied to the ground state 
will crate a superposition of electric and magnetic charges at the neighboring vertices and plaquettes. 

A. Abelian Toric Code 

The simplest choice of G = Z 2 results in Kitaev's canonical toric code [HIES], which has an abelian Z 2 topologically 
ordered ground state. The Z 2 toric code has a 4-fold degenerate ground state on the torus. Since the braiding statistics 
are abelian in this case, the Z 2 toric code can act only as a topologically protected quantum memory and not as a 
universal topological computer [HI UJ3 • 

The Hamiltonian of the Z 2 toric code is given by 

where {cj} are a set of qubits located on the links of a square lattice on a torus, and a" is the a € {x,y,z} 
Pauli operator on qubit j. A v and B p are the 4-body interactions around the vertices v, and plaquettes p of the 
lattice, respectively. The ground state of ifrc is a quantum liquid state in which all local correlation functions decay 
exponentially. Nevertheless, the ground state shows Z 2 topological order, manifested in the four-fold ground state 
degeneracy of H TC on the torus. These degenerate ground states may be distinguished by the action of the following 
non-local loop operators: 

w h = n **> w u = n ° z r ( ie ) 

Here are loops through the vertex lattice that wind around one of the two directions of the torus, and ci,2 are 
loops passing through the faces of the plaquettes. Since [Wf^j-^Tc] = 0, the ground states may also be written as 
eigenstates of W^ 2 Z w ^ n eigenvalues ±1. 

The excited eigenstates with A v \ ^ ex ) = — \^ ex ) have a localized "electric" e-type quasiparticles on the vertex v 
that costs energy 2A e . Eigenstates with B p \^ ex ) = — \^ ex ) have a localized "magnetic" m-type quasiparticles on the 
plaquette p that costs energy 2A m . Pairs of quasiparticles can be created by applying one-body Pauli operators to 
the ground state: 

\e v e v ,) = <T* y |*J C ) , \m p m p ,) = a^ p , \^ c ) (17) 

where cr v y is spin on the link connecting v and v' and cr p p i connects plaquettes p and p' . The lowest energy excited 
states are characterized by such pairs of localized quasiparticles and are separated from the ground state by an energy 
gap of A e , m = 4A e!m . An arbitrary one-body operator acting on an eigenstate of i?TC will act to either create or 
destroy pairs of quasiparticles around the link: sequential action of one-body operators therefore results in translation 
of a single quasiparticle across the link. Both e and m quasiparticles act as bosons under exchange amongst their own 
type: however, braiding an e around an m generates a phase of —1, and so the two types of quasiparticles are seen to 
have mutual abelian semionic braiding statistics. 

On a finite-sized lattice of linear dimension L, the crossover temperature below which the topological order is 
preserved for the Z2 toric code is given by T* ~ A/lnL [UJ. The anyonic braiding statistics allow this abelian 
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toric code to therefore act as the basis for a robust quantum memory as long as it is kept in equilibrium at a low 
temperature, T <C T* . However, as noted above, the abelian topologically protected braiding operations are not 
universal for computation in the Z 2 model. 



B. Nonabelian Toric Code 



Mochon |36l I37| has shown that for certain non-abelian finite groups G, the braiding and fusion of electric and 
magnetic fluxes can lead to formation of a universal logical gate set. Thus, if the corresponding toric code Hamiltonians 
Hq can be robustly generated and localized quasiparticles controllably created and manipulated, the non- Abelian toric 
code may act as the basis for a topologically fault tolerant quantum computer. The smallest non-abelian group is the 
group of permutations of three objects, S3 and this would already allow for a topologically protected universal gate 
set [SHIEZ]- Since S3 has 6 elements, simulating Hs 3 therefore requires using qudits with d = 6 at each link of the 
square lattice. In Section|V]we indicate how the corresponding vertex and plaquette operators of Eqs. ( 10 >-( 1 1 1 may 
be constructed from a universal set of one and two-qudit gates. 



C. Finite Temperature Behavior 



As noted above, while topological order is unstable at any finite temperature in the thermodynamic limit, in a finite 
sized system there is a finite crossover temperature below which the topological order may be preserved [101 SI] • To 
keep the system in such effective low temperature state with respect to Hq, or to cool all the way to the ground state, 
we must therefore couple each vertex and plaquette to a set of ancillary reservoir qudits undergoing dissipation [381139], 
In the next section we describe how this thermalization may be carried out in an efficient manner, focusing on the 
qubit case, i.e., on the abelian Z2 toric code. 



IV. THERMALIZATION OF STABILIZER HAMILTONIANS 



Since noisy control operations and interactions with the environment will introduce entropy to the system, a robust 
quantum simulation requires a dissipative process to remove entropy and effectively cool the system [281 129] . As an 
alternative to algorithmic error correction, we present an approach here to maintain the system in equilibrium with 
an effective thermal reservoir. If the effective temperature of the reservoir approaches zero, this thermalization will 
act purely as passive error correction; however, our approach also provides for the ability to tune the temperature 
of the reservoir and equilibrate the system at finite temperature. In addition to providing robustness, this thermal 
equilibration also therefore would allow for a study of thermal properties of the quantum model at hand. One 
procedure for extracting finite temperature properties of a many-body system is to simulate its thermalization by 
coupling it to external dissipative modes. The steady-state of the following Lindblad master equation is guaranteed 
to produce the thermal equilibrium density matrix of the system: 

f t =~llH a ,p}+J2 ^[Hk]p + E e-^7kV[H+]p (18) 
k k 

where T>[A]p = 2ApA^ — A^Ap — pA\ and H£ creates an excitation of energy Hujk- The unique steady state of this 
master equation is the thermal density matrix corresponding to the inverse temperature /3 and Hamiltonian H |42j . 
Despite this simple description, the mathematical model above can be difficult to simulate in practice because the 
excitation creation and annihilation operators , are usually a superposition of non-local many-body operators 
in an interacting system. For the abelian toric code these are 4-body operators, and therefore each generator in the 
Lindblad terms is a 4-body term. 

These 4-body operators can be implemented using the stroboscopic technique described in Ref . [3U] . However such 
an implementation is resource intensive. Here we describe a simpler implementation of thermalization in stabilizer 
Hamiltonians by exploiting the local properties of excitations. Our analysis in this section focuses on stabilizer 
Hamiltonians of qubits for simplicity, however the arguments presented here may be extended to qudit Hamiltonians, 
as required for realization of, e.g., a non-abelian toric code. The key insight for constructing a simple thermalization 
scheme is to note that a local Pauli operation in the interaction picture defined by the stabilizer Hamiltonian (an 
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eigenoperator decomposition of the local Pauli) decomposes into a small number of Fourier components: 

M 

ei^a^e-i^ = £ e^^a^ + e* 2e — (19) 
k=l 

where a € {x,y, z}, j indexes a qubit in the many-body system, and a a< j t k,a a j k are generalized creation and 
annihilation operators for excitations of energy e a ,j.k, associated with the Pauli operator a". While for a general 
Hamiltonian this eigenoperator decomposition can have a number of terms that is extensive in system size, the 
number of terms above, M, will necessarily be small because the terms in -ff s tab commute with themselves and only 
a small number of them do not commute with the local Pauli operator er". 

In our construction we will interact each local qubit with up to M ancillary spins that are driven to a thermal state. 
The specific form of the interaction Hamiltonian is: 

N M 

^=EEE ff i 8SI ^ ( 2 °) 

j'=l a k = l 

where ^ k is the Pauli x operator on ancilla spin k for the ath Pauli operator on qubit j in the stabilizer Hamiltonian. 
Note that this interaction Hamiltonian is two-body and each lattice site interacts with at most M ancillas.. The 
Zeeman frequency of the fcth ancilla spin is chosen to be u) a ,j,h — £a,j.k/h- 

The effective dynamics of the system and ancillary spins, with combined density matrix £1 e %i attlcc (g> (c 2 ^j® MN : 
is given by: 

-rj: = --[iJ stab - ^2 Hj^^j t + Hi nt ,Q] 

+ E -ra^iPP^kP + E < jt ^K^ (21) 

The dissipation rates for all the ancillary spins are chosen such that 

2^k = ^ptiu^H ( 22 ) 

so that in the absence of the coupling to stabilizer Hamiltonian qubits, the unique steady state of dissipative dynamics 
of a single ancillary spin is the thermal density matrix at inverse temperate /3: 



Pa,j,k 
P**.k 



(23) 



with p\ j k = 1 — p° a j f. = e~P nuJa >i' k p° a ■ k . Such a thermally driven ancilla system can be implemented as a driven 
three-level lambda-configuration atom with a strong spontaneous emission channel, resulting in the two stable levels 
having thermal population distribution; in such a case they are pseudospins. The inverse temperature of the ancilla 
pseudospins is chosen to be j3 — l/fc^T (where T is the desired temperature for the simulation). 

Transforming into the interaction picture with respect to H st&h — ^ Q . fc hu a .j k Y> z a ■ k , and using Eq. ( 19 ) yields an 
interaction picture interaction Hamiltonian of the form: 

N M / M \ 

= E E E E + ^' lt <U ® (S^,,e- 2 — * + K^- ht ) (24) 

j=l a k=l \l=l / 

Remembering that hio a ,j ik = e a .j,k and using the rotating wave approximation to drop all oscillating terms reduces 
the interaction Hamiltonian to 

N M 

^ A = E E E W ® K 3 M + a lj,k ® ^ 3 .k (25) 

j=l a k=l 



Therefore the effective evolution in the interaction picture and under the rotating wave approximation is given by the 
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Figure 1: Graphical representation of the action of excitation creation and translation operators for the toric code. The black 
dot indicates site (qubit) j, the blue (red) circles indicate electric (magnetic) excitations. 



master equation: 



dt 



N 



A I 



^XXXw -" ! 

j = l a k=l 

+ X TaJ^PaJ,*]" + X ^J^Kj,^ 



(26) 



ct,j,k 



a,j,k 



This evolution describes energy exchange of the stabilizer Hamiltonian system with a set of thermalized ancilla spins, 
and will result in the thermalization of the system. Essentially the engineered dissipation channels provided by the 
ancilla pseudospins mimic a bath satisfying detailed balance at the resonance frequencies of the system. Although we 
can demonstrate the thermalization of the system in this general setting (see Appendix) , for concreteness we will now 
illustrate it explicitly for the abelian toric code. 

For the toric code, the Fourier decomposition of the local Pauli terms of interest is: 



x/z - 

' a i e 



-2-| At J7T 

e h E 



3,£/M 



e2lAtE )x/M +T ^£/M 



(27) 



v = £ and v = M. denote electric and magnetic excitations, and we assume the electric and magnetic excitations have 
the same energy: A e = A m = A. The operator ijj v creates a pair of excitations of type v about site j and Tj >v = Tj v 
translates a type-^ excitation about site j. See Ref. [3D] for the formal definition of these operators, and see Fig. jl] 
for graphical representation of the action of these operators. In terms of these operators, the toric code Hamiltonian 
may be written as, 



H 



TC 



X 



2E\ ..Eh 



rp2 



(28) 



This decomposition of toric code excitations was used by Alicki, Fannes, and Horodecki 38 s 39] to show that 
thermalization of the toric code is guaranteed under a local interaction Hamiltonian coupling to a thermal bath of 
the form: 



int, harmonic 



N 

X-i 



N 

X^i®* 



(29) 



where fj,gj are Hermitian bath operators that satisfy detailed balance. For example, in the case of a bath of free 
harmonic modes, fj,gj could be the sum over displacement operators for the modes. We will instead, show that 
by the above construction one can also simulate thermalization by utilizing ancillary pseudospins as the engineered 



dissipative environment. The interaction Hamiltonian we require for the toric code example is 

N 

Him =Y. l® E A,j,£ + "i ® S A ,j,M + °* ® + 0? ® (30) 

Here, the S^/o ^. are Pauli-a; operators on ancillary pseudospins which are implemented as driven three- 
level lambda-configuration atoms with strong spontaneous emission. Each toric code lattice spin has four types 
of ancillary pseudospin coupled to it (£a,£, ^o,s, ^a,m and Xo.jvi)- The characteristic frequency of all A-type 
ancillary pseudospins is wa = A//i, and all 0-type ancillary pseudospins is loq — (these correspond to the 



two eigenfrequencies in the decomposition Eq. (27)). Transforming into the interaction picture with respect to 
Htc ~ J2j v ^ w a^a j u — Ylj ^o^o j vi w hh v £ {E,M}, and employing the rotating wave approximation produces: 

HT A = E E i*> ® S A tj ,u + E L ® S A ,j, v + T jyV ® Eg,^ (31) 

Therefore the effective evolution in the interaction picture and under the rotating wave approximation is given by 
the master equation: 

^ = - JE ^> ® S A J,„ + 4* ® E A,,> + ^> ® S §,^> fi ] 

+ E 7a^Paj,J« + E Ti^Pij-jn + E 7o ^P^, jn + E 7^Pk>]n = ^ (32) 

i.^ i 3 3 

If we choose 7^ = e~^ jA j^ and 7^ = 7^, which results in the A-type ancillary pseudospins thermalizing to a finite 
inverse temperature /3 and the 0-type ancillary pseudospins thermalizing to infinite temperature (/3 = 0), it is shown 
in the Appendix that the unique steady state of this master equation is the thermal state of the combined system at 
inverse temperature j3. That is, 

-I3H TC ^ p -P{- J2j u ^aSaj^-Ej „ Hoi ^oj,v) 

Ces — g = (33) 

and this is the only state that satisfies this property. 

The above analysis for simulating thermalization focused explicitly on qubit stabilizer Hamiltonians. However, 
the same construction follows for qudit stabilizer Hamiltonians because the critical property of local excitations and 



a decomposition analogous to Eq. ( 19 ) also holds for these. In the qudit case however, there are more than three 
Pauli generators in the group of local generators and thus the Fourier decomposition will be more involved. As a 
consequence the number of ancilla spins necessary to simulate the thermal bath will also increase. 

V. PHYSICAL SIMULATION OF STABILIZER SYSTEM 

We shall now consider the physical context in which the aforementioned thermal stabilizer system is to be simulated. 
Though a number of proposals exist for quantum simulation, for specificity we shall restrict our discussion to arrays 
of trapped neutral atoms. In particular, we consider a set of ~ 250 individual 133 Cs atoms trapped at the sites of 
an addressable, simple cubic optical lattice [33]. The orbital degrees of freedom are slow, and can be considered as 
effectively frozen on the time scales relevant to our analysis. We therefore need consider only the internal atomic 
degrees of freedom, from which we select two hyperfine levels (e.g., \F, tuf) = |4, 4) , |3, 3)) to define a 2-level pseudospin 



system. The Hamiltonian, Eq. (15 1, will be implemented in an interaction picture with respect to the atomic energy 
levels. Additionally, we choose auxiliary internal levels to serve as intermediate states to facilitate optical frequency 
Raman transitions for single qubit operations, as well as a highly excited, n ~ 80, Rydberg level necessary for two 
qubit interactions [44j[45]- Since the pseudospins are localized at the sites of a cubic lattice, one can choose to either 
realize the dynamics on a single plane using a surface code |22l 146] or in a three-dimensional cubic array with toroidal 
boundary conditions realized by SWAP operations. In addition to this set of system qubits on which the stabilizer 
Hamiltonian is simulated, we must include a set of ancillary atoms to serve as a thermal reservoir. This ancillary set 
must be strongly dissipative and will be optically pumped to produce the desired thermal state. 

The analysis of the previous section demonstrates that the thermal properties of stabilizer Hamiltonians may be 
studied by coupling to a dissipative bath of two-level systems. However, it is not possibly to directly implement the 
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Liouvillian on the neutral atom system, so we instead take a stroboscopic approach. This approach applies a sequence 
of local operations to generate an evolution which closely approximates the evolution generated by the dynamical 
equations. We shall begin by considering the generic Lindblad master equation, 

p = C P (t) =Y,AHMt) + X>Z>[iW*)- (34) 

3 k 

Here we have used A[Hj]p(t) = —i[Hj, p(t)]/h to represent the adjoint action of Hj. Evolution under this equation 
generates the time evolution operator, 



L(t) = exp(£i) (35) 



The master equation, Eq. ( 34 1 , is impossible to implement continuously in the neutral atom array, so we approximate 



the evolution by a sequence of local operators using a Trotter expansion, 

L(t) = cxp(0(l/A 2 )) ^cxp(^[^]A<)fJexp( 7fc I?[^ fe ]At)j , (36) 

where At = t/N.Hy choosing a sufficiently large N, the error term can be made arbitrarily small. Each term in the 
master equation may then be simulated independently over the short time At. For a generic stabilizer Hamiltonians, 
however, many of these terms will be multibody. Couplings derived from first principles physical interactions, on the 
other hand, are intrinsically two-body. Exceptions to this are rare, and usually derive from an implicit averaging over 
time or intermediate degrees of freedom |47j . In Ref. |30j . we presented a perturbative method based on the Magnus 
expansion for simulating many-body interactions between qubits by application of a sequence of two-body quantum 
gates. Here we present a non-perturbative method for the simulation of arbitrary n-body qubit gates based on a 
construction in Ref. [JHISH]- This method displays significantly lower overhead in terms of total gate count as compared 
to the Magnus expansion approach. While the method is capable of simulating a generic product of Pauli matrices, we 



shall explicitly consider here only the terms necessary for simulation of the abelian toric code Hamiltonian Eq. (15 1. 
The Trotter expansion leaves us to simulate unitary evolution operators of the form Uz(t) = exp(—i^)a z a z a z a z ), with 
4> = XAt. This evolution may be simulated by the sequence of CPHASE and one-qubit gates shown in Fig. |Vj This 
circuit may be readily extended by single qubit operations to achieve any product of four Pauli matrices. For example, 
pre- and post- multiplying by Hadamard transformations produces Ux = exp(—i<fHT x <T x cr x (r x ) = H® 4 UzH® 4 . 



Simulation of the dissipative terms in Eq. ( 34 1 may be effected through the use of atomic states with strong 
spontaneous emission. This dissipation forms the entropy reducing part of the map. As shown in the Appendix, 
the stationary state of the evolution is thermal regardless of the magnitude of the dissipation couplings, 7^, so 
long as they are positive. We may therefore choose them to be arbitrarily large, effectively replacing the operators 
Y\ k exp(T>[Kk]St) by a reset to the thermal state. This could be accomplished in a number of ways. For example, 
correctly tuned coherent driving to a state with strong dissipation can optically pump the atom directly to the thermal 
state. Alternatively, one may measure the ancilla bit in the computational basis and apply a 7r-pulse conditioned on 
the result of the measurement and a classical, Boltzmann- weighted random bit. 




Figure 2: Ancilla-free quantum circuit nonperturbatively implementing action of many body Hamiltonian. The circuit on the 
left utilizes gates available to the neutral atom array described in the text. Vertical lines indicate Rydberg CPHASE gates 
between qubits indicated by black dots and H indicates a Hadamard rotation. Action of Hamiltonian proportional to ^o-jffjffj 
may be implemented by conjugation of the above circuit with Hadamard operations. 



For the non-abelian toric code we need to generate the generalized 4-body interactions in Eq. (10 1. This requires 
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an implementation of the S 3 multiplication table on C 6 , the single qudit Hilbert space. Given an implementation of 
one- and two-qudit gates, we can generate the required 4-body terms by methods similar to those above for the qubit 
case. Brennen et al [50] have demonstrated the existence of universal sets of one and two qudit gates. It remains to 
construct explicit instances of gate sequences. Arbitrary one-qudit gates may be implemented through entirely local 
unitary actions on a six level qudit pseudospin system. In the case of trapped neutral atoms, this is possible with 
methods similar to the single qubit unitaries. For a pair of levels within one six level system and a corresponding pair 
in one of its neighbors, we may construct a Rydberg blockade, analogous to the qubit case. This allows the explicit 
construction of two qudit gates. Details of this will be discussed in a forthcoming paper. 



VI. SUMMARY 



Motivated by the need to incorporate dissipative processes to remove entropy and provide cooling during quantum 
simulations, we have developed a scheme for efficient finite temperature quantum simulation of general stabilizer 
Hamiltonians. These Hamiltonians are typically characterized by non-local or many-body interactions that are hard 
to realize. The well known toric code Hamiltonian of Kitaev, which allows both abelian representations with qubits 
and non-abelian representations with qudits, is taken here as a canonical example of stabilizer Hamiltonian in order to 
demonstrate the approach. Our method relies on coupling of each physical qubit or qudit involved in the Hamiltonian 
simulation to a small number of ancillary pseudospins that are dissipatively driven to reach a specific temperature. By 
using a Fourier decomposition of local Pauli operations on the physical qubits, we show that we can achieve thermaliza- 
tion by employing only two-body couplings between the physical qubits with a small number of ancillary pseudospins. 
This is a significant improvement over our earlier work [30j . which required that the dissipatively driven ancillas be 
coupled to the physical qubits by similar many-body interactions as those contained in the stabilizer Hamiltonians. 
We illustrated the thermalization approach explicitly for the abelian toric code of Kitaev, where two-body interactions 
with two dissipatively driven ancilla pseudospins are all that is required to achieve thermalization of a finite set of 
qubits evolving under the toric code Hamiltonian. This considerably simplifies the physical implementation with 
neutral atoms trapped in optical lattices, for which we also presented an improved approach to quantum simulation 
of the toric code Hamiltonian. The approach can readily be extended to thermalization of quantum simulations with 
general stabilizer Hamiltonians, in particular to the simulation of the non-abelian toric code Hamiltonian, for which 
the smallest qudit dimension is six. Detailed analysis of such non-abelian simulations with trapped neutral atoms will 
be presented elsewhere. 
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Appendix A: Thermal state fixed by evolution 



Here we will show that the thermal state of the abelian toric code is the unique fixed point of the engineered 
dissipative evolution detailed in section IV Consider in particular, the evolution prescribed by Eq. ( 32 ) . The dissipative 



dynamics for the ancillary pseudospins is particularly simple since it is incoherent excitation at rate 7 + and damping 
at rate 7" of each pseudospin independently The only steady-state of this evolution is the mixed state of each 
pseudospin 



f P °o /0 



1 







(37) 



with p^/q = 7a/o/(^a/o ^A/o)" With the choice of rates given in the main text, this is a thermal state at inverse 
temperature j3 for A-type pseudospins and a completely mixed state for the 0-type pseudospins - i.e. Pa/o/^a/o = 



e -/37to A /o t w here p A /o 



1 



Pa/o- 
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Hence the steady state of master equation evolution must have the form: 

n ss = P® P f N ® pf^ (38) 

SApscudospins Sopscudospins 

with additional conditions on /?, the state of the toric code lattice spins. Here, 

i and p = 2 1 

J V u 2 

Substituting this form into the master equation results in: 
dfi ss % 



PA = 1+e Q 0& ! and po = I 2 » ) (39) 



E ® ^ + E L ® S a,> + ® S^.,, p ® pf w ® pf 2W ] (40) 



dt 

since the dissipative terms all evaluate to zero. Now, expanding out the steady-state form for the ancilla pseudospins, 
and requiring that this time derivative be zero results in the following conditions V j, v: 

p A E j ,uP-p 1 A pE j ,u=0 (41) 
PaEI vP -pI P eI u = (42) 
TjMp - P T JM = (43) 

At this point note that the translation operators Tj )V generate a group whose action commutes with the Hamiltonian, 
ifrCi and is ergodic in each energy eigenspace, i.e., any two degenerate energy eigenstates are connected by a product 
of Tj tV (with the exception of the groundspace, which we shall address later). Because the translation operators 
commute with the Hamiltonian and are ergodic on each energy eigenspace, they may be decomposed into irreducible 
representations as 

^>=E P « T ^> P " ( 44 ) 

n 

where P n is the projector onto the n th energy eigenspace. Furthermore, 



The commutation relationship [Ti^,p] = implies then that, 

= {P m T iiV P m ) (P mP P n ) - (PmpPn) (PnTi,uPn) (46) 

Choosing m = n, we see that P m pP m commutes with all elements of the m th irreducible representation of T. By 



Schur's first lemma [5T], this implies that P mP Pm must be proportional to the identity. If we choose m ^ n, Eq. (46) 
and Schur's second lemma [5T] imply that P m pP n = 0. 

For the groundspace, we examine the commutation relations with the string operators. The string operators can 



be r epre sented as exciton pair creations, translations and annihilations. But the "commutation" relations Eq. (41) 



Eq. (43) imply that p commutes with any product of T, E, and for which there are equal numbers of £"s and 



E^s. Because PopPo commutes with all of the string operators, Schur's lemma again implies that it too must be 
proportional to the identity. 



That the populations must satisfy the Boltzman distribution is insured by Eq. (41 ), and so 



Given that the thermal state is the unique steady state of this Lindblad evolution, it can also be shown that it is an 
attractor, meaning that all states converge to it asymptotically [31] . In fact, the above is an explicit demonstration 
of a very general statement about semigroups to be found in the work of Arveson [52] - roughly, if a semigroup 
dynamics (e.g. generated by a Lindblad master equation) has an invariant state, and is ergodic, then it is the 
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unique invariant state and is furthermore an attractor. The ergodicity of the dynamics is the key element, and 
for a stabilizer Hamiltonian it can be shown that if the Lindblad generators are excitation creation, annihilation and 
translation operators the system is ergodic. Such an argument can be used to prove that in the general qudit stabilizer 
Hamiltonian case a construction analogous to the one in Section |IV| will fix the thermal state, and only the thermal 
state, of the system. This reflects a general pattern in the thermalization of stabilizer codes which we will discuss in 
a forthcoming paper. 
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